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ABSTRACT 


Various procedures for the off-line identi- 
fication of a Batch Grainding Mill in noisy environment 
are discussed. A computer simulation of the mill with 
the identified parameters reveals the suitability of 
these schemes in a practical situation. A comparitive 
evaluation is then carried out and conclusions drawn. 



CHAPTER I 


INTRODUCTION 

1.1 Outline of the Problem 

Modelling and identification of parameters of 
the system so modelled is normally the first step in the 
theoretical investigation of the performance of a physical 
system and is of cruicial importance in any optimization 
study or adaptive control of the same. Identification cells 
for input -state or input-output records of the relevant 
system and therefore involves experimental work followed - 
by theoretical analysis. 

This thesis is concerned with a specific problem 
in Metallurgical Engineering which is deemed to belong to 
the area of identification. The system considered pertains 
to the model of a batch grinding ball mill which receives 
dolomite stones of specified size range and as the mill 
revolves, these are ground to finer sizes by the crushing 
and grinding action of steel balls which are mixed with 
the feed. The grinding action is assumed to be described 
by an apriori known set of linear differential equations 
as discussed subsequently. Experimental measurement of the 
distribution of particle sizes for different revolutions of 
the mill represents the data. Although dynamics of the 
plant and the measurement process can only be satisfactorily 
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characterised statistically, as is usual in many a physical 
situation, no statistical information concerning these is 
available. Since the structure or topology of the model and 
the nature of variations of parameters are assumed, the 
problem reduces to one of identifying these parameters to 
satisfy a chosen performance criterion. 

This problem is characterised by the following 

features : 

1. The dynamics of the mill are characterised by 
a continuous model, but measurements are discrete. 

2. Strictly speaking, 2 sets of parameters., the 

so called "selection functions" and the "breakage functions" 
are to be identified. Though the model is linear in state 
variables, the measurements are non linearly related to 
the parameters. This thereby leads to the nonlinear para- 
meter observability problem. 

3. The data available is meagre and ia not regularly 
spaced along the time axis. No ensemble of data is available 
and so the accuracy of any one reading is quite limited. 

To be fair, there are many practical problems in collecting 
accurate data at closer intervals as this is the case where 
measurements disturb the process. 

1.2 Outline of the Dissertation 

Different schemes for the identification of the 
system parameters have been considered in the following 
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chapters. On the whole it has not been possible to give 
any analytical proof for the convergence or uniqueness of 
the procedures employed and numerical convergence has been 
considered satisfactory. 

Ohapter II considers the- mathematical model of 
the grinding process in terms of the selection and breakage 
functions. The model is derived using phenomenological 
considerations and the selection and breakage parameters 
themselves are affected probably in a complicated way by 
many factors like the properties of the material ground, 
the size of steel balls, the proportion of the balls in 
relation to the feed mass and so on. It is not an easy 
task to take into account any of these in our problem. 
Essentially we Judge the success of the modelling by the 
closeness with which the experimental data could be simulated 
analytically. 

Chapter HI covers the first method-Recursive 
least squares. Here the performance measure is the well 
known sum of squares of errors between the simulated model 
and the experimental results. It is sreen that the method 
has not been successful in solving the problem completely. 
While it does give a reasonably satisfactory answer to the 
question of identifying the selection functions, the same is 
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not true as far as breakage functions are concerned. This 
incidentally brings to sharp focus the type of difficulties 
and the lack of any general approach to the nonlinear 
parameter observability problem. What is suitable for one 
problem may not be so good in a slightly different situation. 

Chapter IV discusses the Quasilinearization 
technique as applied to the identification problem using 
the continuous model. The criterion is same as in Chapter 
III and the algorithm provides more or less consistent 
results for both sets of parameters. Finally a complete 
program for the evaluation of the best parameter of the 
system given the' data is presented. 

We discuss the Invariant Imbedding approach to 
the identification problem in Chapter V, The many practical 
difficulties encountered in implementing the estimation 
algorithm for a problem of a specific nature are pointed 
out. Computational results for a set of initial conditions 
are presented and the limitations analysed. 

Chapter VI covers the classical Maximum Likeli- 
hood Identification scheme with Gaussian assumptions. 

Here a more elaborate model with both driving and measure- 
ment noise is identified and is compared with previous results. 
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In addition this method provides an estimate of the various 
noise and error covariances for the assumed model. 

A simple Stochastic Approximation algorithm is 
implemented in chapter VII which is essentially the 
Keifer Wolfowitz scheme with a suitable performance criterion 
The formulation is analogous to the steepest descent technique 
in the deterministic situation. 

Chapter VIII concludes the dissertation with a 
comparitive evaluation of the different methods considered: 
earlier. 



CHAPTER II 


THE MATHEMATICAL MODEL OF THE GRINDING MILL 

Generally speaking, the mathematical modelling 
of the process of crushing and grinding of feed in a mill 
can he extremely involved depending upon the desired 
accuracy and details. Accordingly many models have been 
proposed by various authors (Schuhmann, 1940), (Charles 
1957). In this work the phenomenological model derived by 
conceiving two basic functions in grinding is used. These 
are termed the selection and breakage functions and are 
defined subsequently. This approach has been found to be 
generally satisfactory for characterising and simulating 
of mills of the type considered. In this scheme, a set 
of time differential mass balance equations are presented 
(Reid; 1965) 

dxi(t) i - 1 

~dt = “ S i x i (t) + 2Z % Z 0i x j (t) > 

3 = 1 

i = 1, 2, n (2.1.1) 

"Where Xj_(t) is the mass fraction of particles 
identified by the i th size index at time t and n is the 
number of particle size indices catagorized by two consecutive 
meshes of standard testing screens. S^, the selection function 
is the rate of breakage of particles of i th size expressed 
in mass fraction ground per minute. Zji, the breakage 
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function is the i th size fraction of products when particles 
of size index j are broken. 


in compact form eqn« 
dt 


( 2 . 1 . 2 ) 


where x(t) is the 12 dimensional fractional particle size 
vector and A is 12 x 12 matrix given by 


- Si 0 0 

Sp Zp g — Sg 0 •>••<•••••!• 0 

SpZi 3 S 2 Z 23 -S 3 0 


^ SpZ 112 S 2 Z 212 


“ S 12 


It may be noted that the above equation is a 
convenient approximation to a process which is possibly 
nonlinear and which at the moment does not admit of any 
dynamical noise. An additional error of course is introduced 
at the measurement stage* In this model altogether there 
are 78 unknown parameters - 12 selection functions and 66 
breakage functions. It is felt that this large number of 
unknown parameters precludes any pratloal solution of the 
identification problem without further simplifying assumptions. 



8 


There exists considerable experimental justific- 
ation to assume that the breakage function is normalizable. 
(Herbst & Fuerstenau, 1968), (Brown 1969). This essentially 
means that for particles distributed according to Tylor mesh 
scale it can be assumed that Zjl = ^i - jj i > j. This 
assumption reduces the number of breakage parameters to just 
11 and the total number of unknowns to 23, ' Of the two sets 
of parameters, the selection functions is relatively more 
important, both theoritically and practically. 

The following points motivate a stochastic inter- 
pretation to the problem of identification as applied to 
the grinding mill. 

1. The linear time invariant model discussed 
above is only a convenient approximation to the grinding 
process, and hence is open to question. Accordingly, the 
whole exercise may be viewed as one which fits the best 
model belonging to the class of linear time invariant systems 

fj 

into the experimental data and can be interpreted as a 
parameter optimization problem with an a priori chosen 
performance criterion. 

2. It is not clear whether any additional 
disturbance entering the system dynamically will be able to 
11 explain 11 the data better. The best that can be said here 
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is that it is one possible approach and consequently calls 
for a statistical formulation. 

3. The major factor which justifies a stochastic 
interpretation readily is the presence of possible measure- 
ment errors. Also it is worthwhile to note the readings can be 
critically influenced by speed variations of the mill and 
the exact number of revolutions that has been completed. 

In a practical situation instant acceleration and abrupt 
stopping of the mill are not attainable. 

It is quite possible that most of these errors 
are, to a first degree approximation, zero-mean, and hence 
can be averaged out by an ansemble of measurements by repe- 
ated experimentation.- However, it is noted that such an 
ensemble is not available at the moment for reasons which 
are not entirely obscure. 

The experimental data used in this work for 
identifying the parameters is provided by Berlioz (Berlioz 
L.M. 1966) and is considered to' be reasonably accurate 
among the many data that are available. The particle size 
distribution for three feed masses given by Berlioz are 
tabulated in tables 2.1, 2.2 and 2.3. 



10 


CO 

CO 

H 

CO 

CO 


03 

00 

to 

03 

CO 

rH 

00 

l> 

o 

o 

03 


03 

o 

to 

CO 

to 

00 

O 

03 

rH 

o 

o- 

to 

to 

CO 

CO 

o 

• 

5 

m 

rH 

• 

T — f 

• 

i — 1 

• 

rH 

♦ 

o 

• 

o 

• 

o 

* 

o 

• 

o 

• 



to 

LO 

rH 

CO 

rH 

CO 


03 

o 

LO 

03 

CO 

o 

to 

03 

[> 


CO 

03 

03 

03 

00 

CO 

CO 

rH 

o 

I — 1 

03 

03 

03 

o 

00 

CO 

LO 

CO 

CO 

03 

03 

03 

rH 

rH 

H 

H 

rH 

o 

o 

o 

o 

o 

o 

O 




CO 

03 

CO 

03 

rH 

CO 

to 

o 

03 

o 

rH 

o 

CO 

03 

03 

00 


CO 

03 


rH 


00 

o- 

to 



CO 

j — 1 

03 

o 

LO 


CO 

03 

rH 

rH 

H 

rH 

H 

rH 

rH 

o 

o 

o 

o 

O 

O 

O 

O 


w 

„! 

CO 

CO 

03 

to 

o 

[> 

CO 

00 

CO 


I> 



fell O 

o 

o 

03 

CO 

CO 

op 

o 

03 

o 

to 

03 

H 

3 

o o 

tH 

o 

CO 

03 

CO 



03 

03 

i — ! 

i — { 

H 


M H 

CO 

rH 

rH 

o 

o 

o 

o 

O 

o 

O 

o 

O 

o 


• 

• 

• 

• 

• 

/ # 

• 

• 

• 

« 

m 

« 

o 

t=> 













CO 

hH 













CO 

O 














> 

00 

to 

CO 

03 

CO 

00 

o 

rH 

to 

CO 

l> 

CO 

p\ 

w o 

00 

03 

03 

03 

CO 

03 

rH 

03 

l> 

CO 

03 

03 

w 

00 

00 

o 

03 

00 

to 

CO 

CO 

03 

rH 

! — 1 

O 

o 

m 


CO 

rH 

H 

O 

o 

o 

o 

O 

O 

o 

o 

o 

H 

l-q \ 

• 


• 

« 


• 

• 

• 

♦ 

• 

« 

• 

- • 

hH ! 













03 

H 

Vs-H 1 













1 — I 
kH 



03 

CO 

o 

03 


03 

CO 


CO 

LO 

03 

P3 1 — 1 

o 

CO 

03 

l> 

I — 1 

H 

o 


t> 

H 

03 

o 

co 

♦-H gl 

Jco 

03 

CO 

O 

CO 


CO 

03 

rH 

H 

o 

o 

o 

PQ 



i — j 

rH 

o 

o 

o 

o 

o 

O 

o 

o 

o 

< •• 


' • 

■ • 

• 

* 

• 

» 

• 

• 

6 

• 

• 

• 


CO 

co 

00 

03 

o 

03 

o 

CO 

03 


o 

1 — 1 


03 

03 

00 

03 

CO 

H 

CO 

CO 

uo 

03 

00 

00 


03 

i — 1 

rH 

rH 

O 

o 

o 

co 

' • 

rH 

• 

o 

• 

O 

• 

o 

• 

o 

• 

O 

• 

O 

• 

o 

• 

o 

• 

o 

• 


tn Esjf X 

CD *H 00 

\^m\ 


X 

CD CD 
N T J 
•H £ 
COM 


o 

CO 

03 

CO 



03 


to 

o 

to 

CO 

03 

00 

o 

03 


03 

o 



CO 

03 

03 

o* 

03 


03 

rH 

o 

o 

o 

o 

o 

o 

o 

• 

o 

• 

o 

• 

O 

* 

O 

• 

o 

. • 

o 

# 

o 

• 

o 

* 

o 

• 

o 

• 

o 

• 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

• 

• 

• 

« 

• 

• 

• 

• 

• 

• 

• 

♦ 

rH 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 








o 

o 

o 

o 

o 


^ o 

■a $ 

o ^ 


00 LO 
CO CO 


oo oa 


H CO 

X X 

o a 
o to 
rH H 


H 03 00 ^ tO CD IN 00 03 O H 03 

H rH r-f 



o 

o 

CQ 


11 


03 

03 


H? 

PQ 


ID 

i — I 

to 

LO 

03 

CO 

to 

CO 

00 

00 


LO 

00 

LO 

to 

03 

03 

03 

to 

o 

o 

03 

o 

to 

O 

o 

03 

LO 

00 

i — 1 

CO 

H 

00 

in 

LO 

LO 

O 

• 

o 

♦ 

O 

• 

o 

• 

o 

• 

rH 

• 

H 

• 

H 

* 

o 

ft 

o 

ft 

o 

ft 

o 

ft 

tQ 

to 

to 

00 

00 

03 

CO 

O 

00 

03 

in 

LO 

03 

00 

00 

rH 

o 

o 

o 

to 

o 

00 


03 

00 

LO 

to 

CO 

rH 

H 

o 

00 

to 


CO 

CO 

O 

• 

o 

• 

o 

« 

H 

ft 

H 

* 

rH 

• 

1 — 1 

• 

o 

• 

o 

ft 

o 

ft 

o 

ft 

o 

ft 

o 

o 

in 

T — l 

00 

to 

03 

03 

rH 

03 

00 

IN 

io 

<o 

03 

00 

o 

CO 

00 

to 

03 

in 

t> 

o 

oo 

03 

O 

03 

rH 

03 

IN 

to 



03 

CO 

o 

* 

o 

• 

rH 

« 

rH 

• 

rH 

* 

o 

• 

o 

• 

o 

• 

o 

ft 

O 

« 

o 

ft 

o 

ft 

00 

o 

00 

LO 

rH 

H 

03 


to 


to 

03 

03 

00 

£N 

rH 

CO 


o 

03 

o 

t — I 

to 

00 

03 

oo 

CO 

O 

00 

to 

LO 


CO 

03 

rH 

rH 

03 

1 — I 

ft 

rH 

• 

rH 

* 

o 

• 

o 

* 

o 

« 

o 

ft 

o 

ft 

O 

& 

o 

ft 

o 

ft 

IN- 

o 

i — \ 

o 

00 

H 

00 

to 

cn 

LO 

o 


03 

LO 

Tl* 

o 

03 

03 

L0 

1 — 1 

1 — f 

o 

CO 


03 

LO 

00 

03 

to 

LO 


CO 

03 

03 

H 

rH 

(XI 

• 

1 — 1 

• 

\ — I 

O 

• 

o 

• 

O 

• 

o 

• 

o 

ft 

o 

ft 

O 

• 

o 

ft 

O 

• 

o- 

03 

00 


03 

03 

t> 

00 

IN 


to 

00 

00 

1 — 1 

o 



03 

00 

03 

CO 

CO 

03 

03 

o 

to 

03 

o 

LO 

CO 

CO 

03 

1 — 1 

rH 

o 

rH 

• 

H 

• 

H 

• 

o 

• 

o 

ft 

o 

• 

o 

' ft 

o 

ft 

o 

♦ 

o 

ft 

o 

ft 

O 

« 

o 

oo 


to 

o 

LO 

in 

CO 

LO 

H 

to 

o 

tO 

00 


CO 

LO 

to 

H 

LO 

rH 

03 

to 

00 



03 

LO 

03 

03 

03 

H 

i — 1 

O 

o 

o 

LO 

• 

H 

« 

O 

• 

o 

• 

o 

• 

o 

* 

o 

• 

o 

ft 

o 

ft 

o 

ft 

o 

ft- 

o 

ft 

co 

o 

H 



03 


LO 

CO 

o 


CO 

03 

03 

L0 

03 

00 

03 

o 

o 

LO 

LO 

03 

CO 

03 

o 

LO 

03 

H 

H 

H 

o 

o 

o 

o 

o 

O 

I 9 

H 

• 

o 

♦ 

O 

• 

O 

• 

O 

* 

O 

♦ 

o 

ft 

o 

ft 

o 

ft 

o 

ft 

o 

ft 




o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

ft 

ft 

ft 

ft 

ft 

ft 

• 

ft 

ft 

ft 

ft 

ft 

rH 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 


H OOCO^lOCD D-OOOJO HW 

i — I i — I i — 1 




CO 

♦ 

03 


pq 

H 




O 

i — 1 

to 

03 

L0 

00 

00 

o 

CO 

LO 

to 

o 

CO 



O 

i 10 

1> 


to 

o 

H 


to 

03 

H 

03 

to 



CO 

1 to 

E> 

o 

03 

03 

O 

CO 

to 

to 


CO 

CO 




: o 

» 

O 

9 

I — I 

• 

H 

• 

H 

9 

i — 1 

9 

o 

* 

o 

• 

o 

• 

O 

• 

o 

• 

o 

• 



i 

o 

s 

i> 

03 

03 

o 

o 

H 

03 

LO 

o 

to 

LO 



o 

00 

CO 

00 

03 

o 

to 


to 

03 

03 

CO 



03 

h 

03 

CO 

03 

O 

o 

to 


CO 

03 

H 

H 



i 

h 

1 * 

i 

H 

• 

H 

• 

rH 

• 

H 

• 

o 

9 

o 

• 

o 

• 

o 

• 

O 

o 

• 

o 

9 



o 

i — { 

to 

03 


03 

03 

o 

00 

H 

03 

03 

o 



L0 

LO 

03 

8 

CO 

i — 1 

00 

CO 

CO 

to 

H 


CO 



H 

CO 

to 

H 

00 

LO 

H 

CO 

CO 

CO 

» — 1 

H 



03 

» 

H 

• 

H 

* 

\ — { 

♦ 

o 

9 

o 

9 

o 

9 

o 

9 

o 

• 

O 

9 

o 

9 

O 

• 


CO 















g 















o 














• 

M 


to 

to 

03 

UO 

00 

to 

£> 

to 

LO 

a 

03 

H 

CO 

EH 

O 

o 

03 


to 

to 

H 

03 

03 

t> 

H 

03 

O 

g 

po 

O 

00 

o* 

03 

00 

to 


03 

03 

H 

H 

o 

H 

S 

H 

j — 1 

CO 

1 — f 

H 

o 

o 

o 

O 

o 

o 

O 

o 

O 


o 


0 

9 

• 

• 

41 

• 

9 

• 

• 

* 

• 

« 

o 

> 














to 

pq 














03 

pq 














CO 



LO 

o 

03 

03 

03 

to 

o 

o 

o 

to 

H 

H 




to 

03 

03 

i — l 

CO 

03 

03 

00 


o 

03 

O 

o 

H 

o 

to 

O 

O 

f> 


03 

03 

H 

H 

H 

o 

o 

pq 

H 

00 


H 

H 

o 

o 

o 

O 

o 

O 

O 

o 

o 

pq 

M 

• 

♦ 

• 

• 

♦ 

• 

* 

• 

• 

• 

• 

41 


g 














H 















h 















M 



03 

to 

03 

03 

03 

to 

to 

oo 

00 

LO 

H 

to 

g 



03 

H 

00 

H 

CO 

00 

o 

LO 

o 

CO 

to 

to 



o 

CO 

to 

03 

to 

CO 

CO 

03 

H 

H 

o 

o 

o 

Of 


to 

to 

• 

H 

• 

o 

• 

o 

9 

o 

♦ 

o 

9 

O 

• 

O 

• 

o 

* 

o 

• 

o 

0 

o 

• 


| 














EH 















P 



00 

03 

o 

to 

1 — 1 



H 

03 

00 

to 

03 




H 

CO 

CO 

CO 

to 

o- 

00 

03 

to 

to 

CO 

CO 

to 


o 

O 

o 

o 


03 

H 

H 

O 

o 

o 

o 

o 

o! 



to 

I — 1 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

M 


! 

« 

• 

• 

9 

9 

• 

• 

• 

9 

• 

• 

• 

H 


, 













Pd 















pq 















pq 



o 

03 

o 

00 

CO 

to 

00 

to 

LO 

o 

o 

1 — 1 




to 

00 

o 

CO 

CO 

00 

to 


CO 

CO 

03 

03 



o 

03 

00 


03 

H 

o 

o 

o 

o 

o 

o 

O 



03 

• 

o 

9 

o 

• 

o 

• 

o 

• 

o 

• 

o 

* 

o 

* 

o 

• 

o 

• 

o 

• 

o 

• 



0 0 

iso 

•H £< 

t£M 


O 

* 

H 


O 

O 


o 

o 


o 

o 


03 CO ^ 


o 

o 


L0 


o 

o 


to 


o 

o 


o 

o 


o 

o 


00 03 


o 

o 


o 

H 


o 

o 


H 

t — 1 


o 

o 


03 

H 


12 




CHAPTER III 


RECURSIVE LEAST SQUARES 


3.1 Introduction 

In this Chapter the first method for identifying 
the selection and breakage functions so as to meet a 
chosen criterion is considered. This is the well known 
least square error measure which, historically is one of the 
earliest methods of estimating states and parameters. It 
Is simple, mathematically tractable and conceptually quite 
appealing. The popularity rests partly on the essential 
fact that no statistics about the disturbances affecting 
the system be known. 

The system description chosen is particularly 
simple, ie, no driving noise in the dynamics of the plant 
is assumed and the errors in data are deemed to be solely 
due to measurement noise. 


3.2 Recursive Least Squares via Discrete Model 

As seen in the last chapter the system is des 

cribed by 


x(t) = A x(t) 
and y(t) = x(t) + v(t) 

where x(t) and y(t) are 12 dimensional state and measure- 
ment vectors of particle sizes, v(t) is the measurement 
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noise and A is 12 x 12 matrix given by 


- % 

0 

0 

0 

1 

3 1 Z 1 

-s 2 

o ....... 

0 

; 

• 

• 

S 2 Zi 

“ S 3 

0 

• 

* 


• 

S 1 Z 11 

es* t 

* 2*10 

S 3 Z 9 

• 

-Sl2 


Strictly 

speaking 

the identification of the plant 


matrix A consists of estimating both S and Z vectors 
simultaneously. This develops computationally into a diffi- 
cult problem. Apart from the larger order of simultaneous 
equations .considered, (23 eqns . corresponding to both S and 
Z vectors), the differing sensitivity of these two sets of 
parameters to the output data and the further accenuating of 
the nonlinearity of the parameters-observations relationships 
create difficulties in convergence. This^f partly overcome 
by a step-by-step process of keeping the Z vector at a 
suitable constant value while iterating for S vector and 
vice-versa. For the success of this procedure it is essential 
that there should be monotonic convergence towards the least 
square solution. 

^*3 Recursive Algorithm for S Vector 

The algorithm for updating S given below assumes 
a constant Z for successive iterations. The essential idea 
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is to start with a nominal value of the parameters and 
develop a n updating prodedure which recursively minimises 
the sum of squares of errors over the entire data in a 
sequential manner. 

In the following no probabilistic description of the 
system initial condition is assumed, in accordance with 
the claim that in the Berlioz data of grinding, the feed 
consists of particle size of index one only. Any error 
implied by this assumption is considered negligible. 

The continuous model of the plant in terms of the 
nominal S» and Z vectors is given by 

x°(t) = A° x°(t) (3.3.1) 

y(t) = x° (t ) 

T 

with x° (t 0 ) = (1, 0 0) (3.3.2) 

where x°(t), y(t) and x°(t D ) are of dimension 12. 

Taking the discrete equivalent, 1 and 2 give, 
for a discretizing time interval T, corresponding to 20 
revolutions of the mill. 

A°T 

x°(k*l) = e x°(k), k=l, ...5 (3.3.3) 

x°(k+l) = e x°(k), k=6,7 (3.3.4) 

, 5,0A°T 

and x° (k+1) = e x° (k) , k=8 (3.3.5) 
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y° (k) = x°(k), k=l, 2 9 (3.3.6) 

Let there exist an admissible lower triangular 
matrix A which is stable such that the system given by 

x (k+1 ) = e AT x(k) 

y(k) = x(k) +v(k) (3.3.7) 

with initial condition x(t 0 ) = x°(t 0 ) / where 

y(k), k = 1, 2 9 indeed is the given experimental 

data, such that the sum of squares of errors, 

9 j 

21 v (k) v(k) is minimum with respect to the admissible 
k=l 

class of it matrices. 

It may be noted that in eqn. 6 and 7 the vectors 
x(k) , x(k+l), v(k) and the matrix A are all unknown. This 
beings out the crucial assumption to be made in a noisy 
identification problem where the parameters and hence the 
state trajectory are not known a priori and we are left 
with the noise corrupted measurements of inputs and outputs 
only. Accordingly we express A in the form A = A°+*A with 
the implicit assumption that A 0 is "sufficiently close" 
to A so that x° (k) Is a good approximation to x(k) for all 
k. Analytically, how close should be A 0 to A is a particularly 
difficult question to answer. It evidently depends on the 
signal to noise ratio, and the sensitivity of the system 
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states to changes in parameter values. In a specific prob- 
lem like this, numerical experimentation offers the best 
means of choosing A®. From eqns. 6 and 7, 

yCk+1) = e ^° + *A) T z o (k ) + v(k+1) 

Expanding and retaining anly the first 2 terms, 
A°Tr 1 

y(k+l) = e [i+AAT x°(k) + v(k+l), k=l,...5 (3.3.8) 

Here again our assumption that Afl. is "small" 
comes into. play. In this context, one particularly distur- 
bing thought is that in our problem the time interval T 
has been fixed already by the availability of experimental 
data, making it beyond our control, and moreover T is not 
constant during the. entire length of process. The eventual 
culmination of this is discussed later on. 

It may also be noted that in the derivation k is 
limited to 5 only because of non-uniform step size beyond 
100 revolutions of the mill. Necessary multiplicative 
factors are used in the computation later on to account for 
this . 


Rewriting 8, 

y(k-KL) = e AOT fl+AATl x° (k) + v(k+l) 

pi*. 

= F° [l+AATj x°(k) + v(k+l) (3.3.9) 


k = 1,2, 


5 
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ie, y(k+l) - F°x° (k) = F°aaT x° (k) + v(k+l) 
where 


AA = 


-AS 


°1 


0 


0 

A SpZp -AS2 0 ......... 

£>SpZ 2 A&2^1 -&S3 ......... 


1 ^SSl 2 !! AS 2 Z 10 .AS3Z9 


0 

0 

0 


1 


-&s 12 


le, y (k+1) - F° x° (k) = H(k)AS + v(k+l) 

where H(k)AS = F°AAT x°(k), k=l, ...5 


(3.3.10) 


Here H(k) may be interpreted as the modulation 
matrix (12 x 12) which relates in a linear way the correction 
AS to the L.H.S. which is a function of the measurements 
and the nominal trajectory only. 


Letting r (k“) = y (k+1) - F°x° (k) and renaming 

v(k+l) as v(k), 

r(k) = H(k)AS + v(k) , k = t r 8 (3.3.11) 

which is the linearized equation sought all along in terms 
Of AS. 


Giving equal weightage to errors in all particle 
sizes the performance index to be minimised is 

J = IE ^ T (k) v(k) 
k=2 
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In the actual procedure this minimization of J is 
carried out sequentially as seen later on. However, for a 
single stage of the process, let 

J (k) = v T (k) v(k) . k=2, ... 9 

= [r (k)-H (k)^S k ] T jr (k) -H<k)AS k ] 

Minimising w.r.t. AS k , the necessary condition to 

dJ (k) 

"be satisfied is — = 0 

d&S-^. 

which yields = [H T (k)H(k)[j H T (k)r (k) (3.3.12) 

provided the required inverse exists. This implies that 
H(k) be of full rank. 

Now, considering a multistage process, let r(k+l) 
be the additional Observation" given by 

r (k+1) ' = H(k+1)AS + v(k+l) (3.3.13) 


If we we re to couple equations .* 11 and 13 with 
AS now denoted by AS^^ to indicate that measurements upto 
the k+lst stage have been taken into consideration, we get 


AS k+ i 



lH(k) f f H(k) 1 

- 1 ] 

H (k) 

T 

r (k) 


[H(k+l)j [H(k+l)J 


H (k+1 ) 


r (k+1) 


(3.3.14) 
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It may be noted that the approach necessitates 
inversion of matrices of ever increasing order at every step 
of measurement and no use is made of our previous calcu- 
lation Computationally the implications are quite serious 

we are solving the entire problem all over again, whenever 
a new set of measurements are available. 

To avoid this situation and to process the data 
sequentially we use the Matrix Inversion Lemma (Sage 1968). 

The lemma is given below and is easily proved. 

If P k i x = P^ 1 + H r (k+1) H (k+1) , then 

p k+l = p k - p k H T (k+l) [ I+H (k+1 ) P k H T (k+1 )] - 1 

H (k+1 ) P k (3.3.15) 

where all required inverses exist. 

A few words about the notation used may be 
appropriate at this stage. When the index k is used as 
an argument, for example H(k), it indicates the variable at 
the k th stage only. On the other hand when used as a 
subscript, for example Pk or ASfc we imply the variable to 
be the result of calculations utilizing all the stages I, 

1 ^ i ^ k. 

Rewriting eqns. 12 and 14, 

AS^. = P k H T (k) r(k) 


(3.3.16) 
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and AS k+1 = P k+1 £ H T (k ) r (k )+H T (k+1 ) r (k+1 ) ] (3.3.17) 


where 


P k X = H T (k)H(k) 

-1 


a nd P k+1 = H (k)H(k)+H T (k+l)H(k+l) 


Substituting for H i (k)r(k) in terms of P k+ i 
obtained by definition of the latter into eqn. 17 we get 


i 


^Sk+l = p k+l [ p k+l'* H ' I ’^ k+1 ^ H ^ +1 ^ AS k +P k4 ql3 T (k+l)r (k+1.) 
e, AS k+1 = AS k +P k+1 H T (k+l) [ r (k+1 ) -H (k+1 ) & S k J (3.3.18) 


Also by using the Matrix Inversion Lemma, 
p k+l = p k“ p k HT Ck+1 ) [i+H (k+1 ) Pj^H T (k+1 )] - 1 H (k+1 ) P k (3.3.19) 

The algorithm is implemented with starting values 
for S and Z andAS; calculated sequentially taking into 
account the entire data. The starting values are now up- 
dated by S <r* S+ AS and the entire procedure rep.aated. 

In every case the nominal solution using the discrete model 
and initial conditions was obtained. The sum of squares of 
errors involved is also calculated and in fact gives the 
clue to convergence or otherwise of the problem. The results 
are tabulated in table 3.1. It is seen that we get fairly 
satisfactory convergence of all selection functions for diff- 
erent starting values. 
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3 .4 Res u lts 

The convergence of the selection functions 
3 1> s 2> Sn and S»3_2 which were considered representative, 
are shown in fig. 3.1 and 3.2 for different starting values. 
The latter also shows the minimization of the sum of squares 
of error for a close starting value, (case I of table 3.1) 

It may he noted that even for this close initial guess, the 
convergence is slow and oscillatory. 

The same approach is next tried for updating the 
value of Z, the breakage function. It is found that the 
procedure is not successful as there is no convergence 
even with close starting values. The behaviour is oscillatory 
with eventual divergence. 

It is apparent that on the whole the data is not 
sufficient to draw any conclusion straightaway and such an 
attemptmay be presumptuous . But numerically the trouble is 
mainly about the lower elements of the Z vector which get 
modified disproportionally resulting in build up of errors. 

It appeals that for the discrete model the sensitivity of 
the different components of the Z vector to the data is widely 
different. It may also be noted that the large and non- 
uniform discretization of the system can introduce appreciable 
error in the algorithm which neglects the higher order terms 
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in eqn. 3.3.8* This is indirectly verified in the quasi- 
linearization technique where we use the continuous model 
and sampling is done at closer intervals. 

To summarise the positive and negative aspects of 
the above method, we have a workable algorithm for the identi- 
fication of the selection functions when there exists a reasona- 
bly good knowledge of the breakage parameters. In this 
connection it is pointed out that the latter can be estimated 
very simply with due assumptions on the precision of the 
readings at 20 revolutions of the mill as follows! 

By definition is the i th size fraction of 

products when particles of size index j are broken. Putting 
3 = 1 and i = 2, 3 ..... 12, and considering a short interval 
of the mill with feed of only size index 1, we get 

x i+l 

Z± - — i * 1,2,3 11, 

1-xi ’ 

both x-j+i and xp measured at 20 revolutions of the mill. 

The above relation generally suffers from the limited 
accuracy of a particular set of measurements and the assumption 
that secondary breakage of daughter particles can be neglected 
in the interval. 
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In fact the readings for the convergence of 3 
shown in table 3.1 correspond to this particular value of 
Z. When compared with the refined parameter estimates given 
by the C?uasilinearization method, what is obtained here 
by the above, it is seen, is only marginally inferior. 
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P co 
O £ 
P 

p 

I — I d I — I 11 

H W 
M (XI 
■PJH 

CD 
s CO 

: ctf 
‘O 

P 

P 

l ^ 

1 -P 


o 

0) 

H 


LQ 

rH 

o 

CO 

CO 

LO 

CO 

to 

CO 

H 

to 

(XI 

O 

l — 1 

o 


T — 1 

CO 

<2> 

o> 

G> 

a> 

CO 




uo 

Oi 

o 

o 

o> 


EN 

CO 

A 

CO 

co 


CO 

CO 

H 

H 

1 — 1 

o 

o 

O 

O 

t> 

♦ 

• 

• 

♦ 

• 

♦ 

• 

• 

• 

• 

• 

* 

rH 


OOOOOOOOOOOOOOCOOOOOCOOO 


o o o 


o o 


O O O CO 

1> 


■ Cm CO 
O £ 
p 

«r5 -P 
d m 

o 

P H 
M <3 


H 

(X) 

o 


a> 

CO 

to 

tH 

o 

o 

rH 

^ 1 

a 

H 

to 

a) 

G> 



00 

LO 

G> 

00 

O 

CO! 

o 

CO 




LO 

<J) 

o 

o 


LO 

00 

CO 

• 

to 

LO 


CO 

(X 

H 

1 — 1 

rH 

O 

o 

o 

°l 

00 

• 

• 

« 

« 

« 

• 

♦ 

• 

* 

• 

♦ 

• 

1 — 1 


00 CO OO’OO CO CO 00 CD CO CO. CO 00 


o o 


o o o 


T5 fn 

00 


to 

CO 

00 

00 

to 

00 

o 

(X 

rH 

CO 

CO 

d -p 

O 

LO 

rH 

to 

rH 

o 


o 

LO 

no 

to 


LO 

H M 

(XI 




LO 


o 

o 

G> 



(X 

• 

to 

LO 


CO 

CO 

H 

H 

rH 

O 

O 

O 

o 

00 

4- > 00 
<i! 

• 

• 

• 

• 

* 


* 

a 

• 

• 

♦ 

• 

H 


LO 

LO 

o 

LO 

o 

LO 

o 

LO 

(X 

o 

00 

00 

o 


LO 

LO 


CO 

CO 

CO 

03 

H 

rH 

rH 

o 

o 

to 

p 

! 

• 

• 

• 

• 

• 

• 

• 

« 

• 

• 

• 

• 

p 

05 

p 

M 

i O 

t 

1 

1 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

00 

rH 


a> a 

1 H 

LO 

0> 

CO 

CO 

<o 


CO 

(X 

LO 



ho o 

LO 

(X 

LO 

!> 

o 

CO 

O) 


to 

CO 



ctf -H 

CO 


03 

to 


LO 

0> 

o 

00 

rH 

o 

1 

M +3 

<x 

H 

o 

to 


co 

! — 1 

03 

rH 

I — t 

! — 1 

1 

CO o 

t 

03 

H 

O 

o 

o 

o 

o 

O 

o 

o 

» 


© £ 

• 

» 

• 

• 

• 

» . 

• 

• 

* 

-• 



HCXOO^tOtOOCO 


H CO 
rH H 


MILL FEED 3300 GMS. 





1 9H 


Emm 




Ililll 


CONVERGENCE OF S 

__ — — — — — — — * — * — — 

FOR DIFFEREN T INITIAL 
MILL FEED 330C 


ERROR X 


CHAPTER IV 


QUASILINEARIZATION VIA THE GONTINOUS MODEL 
4.1 Introduction 

The parameter identification problem except in 
trivial cases being essentially a nonlinear estimation 
problem which even in the determiniatic cases defies straight 
forward analytical solution, it is quite natural to attempt 
linearization about a chosen nominal parameter value in an 
endeavour to successively update the same. This approach 
has been found to be quite powerful in our problem and gives 
more or less very satisfactory results. 

In a broad framework, the quasilinearization technique 
is essentially a generalised Newton- Raphs on method for 
functional equations. The original nonlinear system is solved 
by solving a sequence of linear systems which when conver- 
gence occurs finally gives, the solution for the original 
problem. The motivation for such an approach is that a non- 
linear multipoint boundary value problem is very difficult 
to solve analytically, and is not particularly well suited 
for computer solution. However this is not the case with linear 
multipoint boundary value problems where the powerful super- 
position principle holds good and the problem lends it self 
for effectiv 
algorithm is 


computational methods . The Quasilinearization 
Characterised by two essential features - 
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quadratic convergence and monotonic convergence. Consistency 
of the results is numerically verified by giving widely 
different starting values for the parameters to be identified. 

The application of the method of quasilinearization 
to an identification problem consists of augmenting the states 
by expressing the (time invariant) parameters in the form of 
a trivially dynamic system. The "state equation" is now 
clearly non linear which is expanded in Taylor series about 
the nominal parameter values, retaining only the linear terms,. 

At this stage what we obtain is a time-varying linear system 
with unknown initial conditions. In fact the missing initial 
conditions are determined to satisfy an apriori criterion that 
the solution should satisfy. And in this case the performance 
index is again the sum of squares of the errors over the 
entire data. The problem is formulated using the continuous 
model keeping the discretization interval free to be chosen 
in the solution. In common with the method of least squares, 
we assume no driving disturbance and are satisfied incorporating 
measurement noise only. 

4.2 Derivation of the Quasilinearization Algorithm 

As mentioned in section 3-2, the identification 
of parameters is done in 2 stages which follow each other 
in an alternative fashion. 
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a. An iterative procedure of estimating S vector 
with Z maintained constant throughout at a nominal value to 
begin with, which is subsequently updated by step b. 

b. An almost identical procedure as above, with 

the performance index minimised with respect to Z in a sequential 
manner, keeping S at the best value obtained by step a* 

Algorithm for Step a 

We have the system 

x(t) = A(S) x(t) where the 12x12 matrix A is known in • 
terms of the unknown parameter S, and x(t 0 ) is known* 

y(t) = x(t) + v(t) 

where y(t) is the measurement vector 
and v (t) measurement noise. 


Consider the N + 1st iteration of the algorithm, 
the iteration index being indicated by the right superscript. 

,N+1 n+1 


. IT+1 
x (t) 


ACS ^ x X (t) , N =1, 2, 


(4.2*1) 


Expanding the R.H.S. about the nominal value 
and x N (t), 
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•N+l TvT AT 

x (t) = A(S N )x N (t) + 

d 

( A(S)x(t)T ■ 

[x N+1 (t)-^(t)] 


c)>c 

1 Jx=x N 



- 




{ A(S)x(t) } 

OS 1 J 


N 

x=x 

$=S N 


0 


S H+1 . a*' 


(4.2.2) 


= A(S N )x N (t) + A(S N ) !_x N+1 (t) - x N (t)' 


L [x :i (t), S N ] [s 


N+l _N "I 

— O 


(4.2.3) 


where L jx N (t) ,S N ] = — j~A(S)x(t) 


_ x(t)=x N (t) 


3 =S 


N 


(4.2.4) 


In lieu of the a priori assumption that 3 is 


- N+l 

time independent, S (t) = 0 


(4.2.5) 


Augmenting the states we get 


r.N+i 

X 



.N+l 

s (t) 


(24x1) 


A(S N ) • L(3 N ,t) 


' x N+1 (t)'“ 


L(S N ,t)S N 

■o' - V 


_N+1 . , 

S (t) 


0 

L 

(24x24) 


(24x1) 


(24x1) 


(4.2.6) 



30 


In principle, the R.H.S. is fully known-both the 
parameters and states correspond to the N th iteration; and 

equation 6 can be solved for any initial condition. Formally 
writing down the solution, 




luCt, t 0 ) ; 0 12 ct, t 0 ) 

3 W+1 (t) ! 


r 

0 2 i(t, t Q ) j 0 22 (t, t G ) 


j 


x M+1 (t 0 f 

JL 

Pi(t) 

s ff+1 (t 0 ) 

r 

p 2 (t) 


(4.2.7) 


where Pl (t) and p 2 (t) are particular solutions of 6 and 

2(t, t 0 ) which is shown in partitioned form is the state 
transition matrix. 


As far as the solution is concerned, it is clear 
that p 2 (t) = 0 and since S N+ l(t) = S N+1 (t 0 ) for all time t, 

^21 (t? t D ) - 0 and 0 22 (t, t 0 ) = I. So writing down the non- 
trivial part of eqn. 7, we have 

N+l 

* ft) = 0 u (t, t 0 )x N+1 (t 0 ) + 0l 3 (t,t o )S w+ l(t o ) +Pl (t) 

(4.2.8) 


011 (t, t D ) and 01 2 (t,t o ) 
overall state transition matrix $(t 


^(t, t Q ) = 


A(S n ) 

0 


I 

t 

t 

t 

l 


L(t 


are obtained from the 
t c ) determined by 

$(t>*0 (4.2.9) 


WvtK 


- i_ 
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And finally P]_(t) is given by the following 
differential equation 

Pl(t) = A(S N ) Pl (t) - L(t, S N )S F (4.2.10) 

with initial condition Pi(t 0 ) = 0 

Now the solution takes the form 

x 1 (t) = 0u(t,t o )x W+1 Ct o )+0 l2 (t,t o )S N+:L Ct o )+p 1 Ct) (4.2.11) 

At this stage, the essentials of what follows can 
be summed up to see where all this ultimately leads. It 
may be emphasised here that Sr I+1 (t 0 ) is still unknown. 

In principle 0^i(t,t o ) and 0i2(t,t o ) and p^(t) 
can be solved from the known initial conditions and x N+1 (t) 

for any t can be obtained in terms of the missing initial 
conditions S N+1 (t 0 ) # 

Let these equations be of form 

x N 1 (t 1 ) = H-j_ S N+ 1 (t 0 ) +w(t 1 ) 

• • 

• « 

♦ • 

x N+ 1 (tg) = Hg S I,+ 1 (t 0 ) +w(t 8 ) (4.2.12) 

where t 1? t 2 ...... t 8 correspond to the time instants at- 
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which measurements are available and H^j Hg ..... Hg are 

completely known. Hence assuming measurement errors, the 
problem reduces to the determination of Si^ + ^(t 0 ) so that 
x^ + l(tp), ........ , x^' + ^(tg) fit into the given data 

in the least squares fashion. 

It can be easily seen that the forgoing exercise 
will be computationally very tedious and time consuming. 

To briefly summarise these requirements we have to solve 
two sets of 12x12 simultaneous (one of them time varying) 
linear first order differential equations for the state 
transition matrix, and another set of 12 linear time- 
varying first order differential equations for the parti- 
cular solution p 1 (t). This will enable one to determine 

Hp, i = 1, 2 8 to further proceed with the least 

square evaluation of S^ +1 . 

It was not considered desirable to proceed as 
above, but instead exploit the condition that 3^+1 (t) is 
time invariant and can be related to S^(t) by a simple 
constant “error term” AS, which is consequently evaluated 
using essentially all the forgoing ideas. The following 
accordingly eliminates the necessity of calculating a time 
varying state transition matrix, but employs a numerical 
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integration procedure which is conceptually and computa- 
tionally straight forward. 


-N+l 


we have x"' 1 ' = A(S N+1 ) x N+1 (t) 


A(S> N )x N (t) 


— (A(S5)x(t) j 


x=x 

3=3 


x N+1 (t) -x N (t) 


as 


y 

A(S)x(t) 


x=x N 
S=SN • 


I N+1 - S N 


= A(s N )x N ct) + ACS N ) 



(t) - x N (t) 


+ L(x N (t), S N ) AS N (4.2.13) 

where L(x N (t), 3 N ) is given by eqn. (4.2.4) 
ie ? 1 (t) = A (S N )x N+1 (t) + L(x W , S N ,t)^S^ (4.2.14) 

with x N+1 (t 0 ) = (1, 0 .... 0) T . 


Let $(t,t Q ) be the state transition matrix for 
the time invariant system 14 given by 

S(t,t 0 ) = f(t-t 0 ) - /(S N )(t-t 0 ) (4.2.15) 

with $(t 0 ,t 0 ) = I 
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1 

x N+1 (t) = $(t-t 0 )x N+1 (t 0 ) + J )L(f) A&: dt (4.2.16) 

„ t Q 

At t = t]_ wVx-QY-e L. <s v* (f cNi > 

+ ti 

x W+ ' 1 (t 1 ) = $(t 1 -t 0 )x N+1 (t 0 ) + I 5(t 1 -T)L(t) Asf'dt* 


le, x N+1 (t 1 ) - $Ct 1 -t 0 )x iA1+1 (t 0 ) = R ± AS 


N+l , 


(4.2.17) 


Considering the observations y(t) = x N+1 (t)+v(t) 


we rewrite 17, 


y(t 1 )-v(t 1 ) - 5(t 1 -t 0 )x N+1 (t 0 ) = H 1 a>3 


Minimising 


v(t]_) 


R 


-1 


where R is the eovariance 


of the observation error gives, 





T -1 


^ y(ti)-$(t 1 -t o )x W+1 (t 0 ) 


(4.2.18) 


Considering subsequent observations at t k , k=2, ....8, 
we can obtain AS^ sequentially by the use of the matrix 
inversion lemma as in section 3.3. The sequential estimation 
equations are i 


P k+1 = P k~ P k H T (k+l) 


R(k+l)+H(k+l)P k H T (k+l)l“ 1 


J 


(4.2.19) 


H(k+l)P k 
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ASS k+l = + ' p k+i H x (k+1 )R“ 1 (k+1 ) 


y ( \+i } (tk +1 -t 0 ) X N+1 c t 0 > 


- H (k+1 ) A3 k 


(4.2.20) 


Algorithm for Step b. 

Here the S; parameters are kept constant and updating 
of Z parameters is done. The essential ideas remain as in 
the previous algorithm. 


x N+1 (t) = A (Z N+1 ) x N+1 (t) 


- A(Z^)x^(t)+ 


lr{ A(z)x(t } - 




‘ f 




A(Z)x(t ) 


_ N 

x-x 

Z=Z N J 


* x=x N 
Z=Z 1 ' r ~ 1 


N+l N 

z - z 


x W+1 (t)-x N (t) 


.N+l 


ie, x (t) = A(Z^)x^ + ^(t) + M-o. Z 


(4.2.21) 


where M(Z N , x N (t)) 


a 


A(Z)x(t) 


x=x 


N 


Z=Z 


N 


The procedure is exactly similar to what was 

done previously. The H(i) matrices i = 1, 8 are 

of dimension Hxll. 
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4.3 Computation Considerations 

Each iteration consists of the following steps. 

1. Nominal solution of the system with assumed 
or updated values of the parameters. The discretization 
interval was chosen quite small (corresponding to 2.5 
revolutions of the mill) to evaluate the integral in eqn. 
4,2.16 numerically. 

2. Determination of H(i), i = 1, 2 ...... 8 

and sequential estimation of i>, (or £>. as the case may 
be)using equns. 4.2.19 and 4.2.20 (or corresponding eqns,. ), 

3. Updating the starting parameter values and 
return to step 1. 

As was indicated earlier, the entire problem 
was tackled by alternating steps a and b. For this approach 
to be successful, it is necessary that we get monotonic 
convergence in both the steps a and b. This is first of all 
verified for different starting values for both S and Z 
parameters to be true. Consistent values for the parameters 
and the sum of error squares -toere obtained for widely different 
initial guesses. These are tabulated in tables 4.1 and 4.2« 
The behaviour of S 1? Sg, Sqq, S 12 , Zq, Zg and Zqq are shown 
in figs. 4.1, 4.2 and 4.3. 
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The same program was used for alternate iteration 
of S and Z vectors with suitable control statements and 
hence a complete pack for the identification of all para- 
meters in the least square fashion was developed. By a 
study of the rate of convergence of the algorithm it was 

decided to have 4 iterations per step for both parts a and 

s^tis |-acfcory 

b, which alternate each other. This was considered of^rrmra i 
to attempt overall rapid convergence. It was also arranged 
such that the program automatically switches itself from one 
step to the other if the '’gain' 1 in the performance index is 
smaller than a prescribed figure at the end of every iteration. 
The stopping criterion is a gain of less than 1 % in the 
performance index for both 3 and Z parameters, or 16 itera- 
tions, which takes almost 8 minutes in Fortran. If necessary 
the program can be rerun with the last updated values of 
parameters as inputs. But this will be unlikely if the 
initial guesses are not off by more than 20$. 

It is also desirable from the point of view of 

np. sfr 

lesser^iterations that the initial starting value for Z is 
calculated by the approximate procedure given in chapter III. 
The results for the different feed masses are given in 
tables 4.3, 4.4 and 4.5. 
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Considering the accuracy with which the data 
has been fitted, we get the sum of squares of error for 
3300 gms feed to be .001428 whereas the sum of squares of 
the data itself is 2.16401, representing an average error 
of 2.6%. This seems to be a considerable improvement over 
the values which were known previously. (Average error 
5.9%). The procedure was tried for feed masses 1980 gms. 
and 3960 gms. with resulting average error 2.9% and 3.4% 
respectively. The mono tonic convergence is shown for all 
the cases in fig. 4. 4. 

4.4 Discussion 

We considered the quasilinearization method 
for solving a multiple point boundary value problem in the 
least square manner. There are many points in this chapter 
in common with chapter III. The criterion is the same, 
and we essentially try to force the nominal solution fur- 
nished by the assumed parameters to approach the data 
in a_ least square manner by a recursive procedure. 

Even the first method implies a sort of linearization by 
expanding e^ T and retaining only the linear terms. The 
essential difference between the two methods is the nature 
of the system eqns. tackled. In the first we formulate 
the linearized system in discrete form^ where the time 
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interval was beyond our choice. And mostly due to this, 
the first algorithm is not successful for the iteration of 
Z vector. This flaw is removed in the Quasilinearization 
algorithm which has rapid and monotonic convergence for 
oth & and Z parameters. As far as the results are concerned 
it is noted that quasilinearization is decidedly superior 
in as much as convergence is monotonic. We get consistent 
results which compare favourably with that of chapter III. 
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CHAPTER - V 

?iB ^-g^S°Biirop3 imR-rawr dibedding 

5,1 Introduction 

The mathematical model describing the entire 

Process of grinding chosen in the previous chapters was 
Particularly slmple . lt dw mt ^ of any drivlng 

noise affecting the syst m in a dynamical way . u lg true 
that so long a s we do not have any knowledge of noise 
statistics, it is Material to put them at the driving end 
or measurement end of a linear dynamic system. Nevertheless 
it is interesting to conceive of a disturbance input for 
the system and reformulate the problem in a different way.,. 
This is done utilising the Invariant Imbedding approach 
for the identification problem as discussed below. 


In optimization studies, we generally encounter 
a two-point or multi -point boundary value problem where the 
terminal conditions are split. It may be emphasised here 
that the concept of invariant imbedding is quite general 
and powerful. Whereas quasilinearization represents an 
iterative procedure where a sequence of linear problems 
are solved, invariant imbedding tackles the original 
problem by expanding it into a family of problems . For 
example instead of considering a process of duration T 
the invariant imbedding formulation is to consider a set 



*±o 


of problems 
zero to T, 


with the duration of the process 
These problems are then imbedded 


particular original problem. 


ranging from 
to obtain the 


fication is 


The invariant imbedding algorithm for identi- 
foimulated in terms of the continuous grinding 


model of chapter II. 


As can be seen later, in the algorithm 


that we obtain, 


by virtue of the moderately large dimension 


the system variables and unknown system parameters, 
the solution is very cumbersome in terms of the number of 


imultaneous non-linear differential equations to be solved., 
encs m order to limit the computational requirements to 


reasonable levels it is assumed that the breakage functions 
are known within close range and identification of selection 
functions only is called for. 


One more aspect of the problem, ie, the non- 
availability of readings along the time axis at uniform 

intervals is of great relevance here. In order that the 

3 

algorithm may be implemented with reasonable accuracy, it is 
imperative to have measurements at closer intervals . than 
is available in Berlioz data. Hence it becomes alm ost 
unavoidable to resort some from of interpolation of the 
given readings. 
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^ ^ ThQ I nvariant Imbeddin g Algorithm 

As in the quasilinearization technique we 
augment the states of the dynamic system by expressing the 
constant selection parameters as trivially dynamic. In 
the following derivation the vector x represents this augmented 
state of the system. Hence we now have 

&(t) = g(t, x) + k(t, x) u (t) 

y(t) = h(t, x) + v(t) , 0 ^.t (5.2.1) 

where x is a n vector 

y is m vector of observations 
v is m vector of measurement error 
u is unknown noise assumed to be a r vector 
g(t,x) is a known n vector function 
h(t,x) is known m vector function 
k(t,x) is known nxr matrix function 

To obtain a logical performance index we define 
the vector valued residual errors as 

e^Ct) = y(t) - h(t, x) (5.2.2) 

62 (t ) = x(t) g(t, x) (5.2.3) 

where x(t) is the estimate of x(t) 9 O^t^T. 
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Using the least square philosophy, the per- 
formance index is now 


J 


J 

0 


e 1 (t) 


+ II e 2 (t) 


dt (5.2.4) 


W 


which is to be minimised with respect to x(t), 0 St . 

Q(t, x) and W(t, x) are at least positive semidefinite 
weighting matrices. 

Substituting equations 1,2,3 into eqn. 4, 

T 


J — 


0 


■ 2 

y(t) - h(t, x) | + |u(t) | 

*9 ' V 


dt 


, T 
k Wk 


(5.2.5) 


where the minimisation is now with respect to both x(t) 
and u(t), O^-t^T subject to the differential constraint 


x = g(t, x) + k(t, x) u(t) 

Let v(t, x) = k T (t, x) W(t, x) k(t, x) 

Define the Hamiltonian H(t, x, p, u) 


(5.2.6) 


by H(t,x,p,u) = |y(t)-h(t,x)J| 2 + 


* V 


\P? g (t,x)+k (t , x}u^> 
(5.2.7) 
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-here p is the n dimensional oo-state variable. 

The unknown input S can be eliminated from 
eqn. 5.2.7 by the Pontryagin's Minimum Principle 


b H 


<3 U 


0 


which should be s&ti of** 

e satisfied along the optimal trajectory 

which is now denoted x* and the 

as H*. 


corresponding Hamiltonian 


Carrying out the elimination we get 

H*(t,x*,p») =|y(t)-h(t,x*)| Q 8 + ^, g (t,x*>) _x <p*, kl rVp) 

(5.2.8) 

Now the reduced state costate equations that 
the optimal trajectory should satisfy/the following canonical 


equations 


. * H* 

x * ~ Tv*~ (t? x *> p¥) 


p* = 


~aH* 

d x* 


(t,X*,p*) 


(5.2.9) 

(5.2.10) 


Since T has been fixed an x*(0) and x*(T) are free, 
the transversality condition yields 

P*(0) = 0 , p*(T) = 0 


(5.2.11) 
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reduces to ^ SeqUentlal estlmatl °n P^obl, 

VI ng the two point boundary value 

P T Chara0tenSed * -- 9 - “ - boundary ' 

™ onsofll> It may be noted that both x*(0) and 

^ arS Unkn0 ™- 1116 ^^ing approach ls 

1Sld8r a largSr ° laSS 0f TPBVP ' S with boundary 
conaitions # 


P*(°) - 0 and p* (^f) 


(5.2.12) 


where T is recognized as the running time variable and c 
n vector. Now the missing terminal condition on x ls 
denoted by r(c,T) where r is unknown. It can be easily be 
s own by perturbation techniques that r(c,J) satisfies the 
nvariant imbedding equation (Detehmandy 1965, Sage 1968) 


£ r ^ H* 

^ G > r 


(T. r* o') - ^ H* 

’ ) c : (T > r ) °) (5.2.13) 


where J-— -I 

I <*c I. 


Hi 

^ Ci 


The above nonlinear partial differential 

equation is approximately solved by assuming a solution of 
form 

r(t, T) = p(i ) 0 + S (T) 

(5.2.14) 
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where P(T) is an n x n matrix 

Substituting eqn. 14 into eqn,, 13 , 

dp 


• dT 


c + 


d x 


- P(T) - H * t m p - . d H* 

T dr CT ’ Pc+ x 5 c)=: •— (T, Pc+x, c) 

(5.2.15) 


sxtj wq expand eqn. 15 about p (o t 1 ! r>a-t- 

ho aDour rw,TJ retaining terms 

order. The justification for this approach is that 

only those solutions of eqn. 13 f or which c=0 are of lnt6rest 

to us. Ana the least square estimate of x(f) denoted by id) 

1$ r(0 ’ T) - After ending eqn. 15, we collect terms of 
order e»(ie, independent of c), and e. The desired sequential 
estimator equations •can be found by equaling coefficient of 
c° and e to sere because c is completely afbitrary. After 
nsiderable algebra we get (Detchmandy 1965) 


x(T) = g (x,T)+2P(T) 


dh^(x,T) 


x 


Q [y(T)-h(x,T) ] (5.2.16) 


?(T) 


^ g(x, T) ^ e T (x T') 

— — P(T) + P(T) — ~ 


+ 2P(T) 



dh T (x,iy 

Q /y(T)-h(x,T)) 

X ^ J 


+ i k(T, x) Y -1 (T, &) k T (T, x) 


P(T) 


(5.2.17) 


1.1. T. KANPUR 

CENTRAL LIBRARY 


.211 PS 
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Equations 16 and 17 represent the algorithm 

f °V° Se9Uentlal eStl “ tion the states and parameters 

° e . SySt : H ’ taklDg “count measurements one by one. 

o major difficulty inherent in this approach is that the 

initial conditions x( 0 ) and p/m 

ana P(0) are not known and is 

subsequently disau^^/nri ,, 

y uascussed m the next section. 

5 * 3 £^E j^ tion al Considerations 

^ Slnce thQ eont “uous imbedding formulation of 

the problem is utilised, we require data that is more 

closely spaced than Is available in Berlioz -s. A linear 
interpolation of the data is affected to obtain readings 
at every 5 revolutions of the mill. 

ing to the implementation of equns . 16 and 
17, what we have is a pair of coupled multidimensional 
nonlinear ordinary first order differential equations. 

For a n dimensional system with m unknown parameters, a 
total of (n..+ m) + (n + m) 2 coupled nonlinear first order 
differential equations have to be solved with unknown 
initial conditions. Hence computationally this method is 
not well suited for a problem of even moderate dimension. 
Convergence is not assured for an arbitrary starting guess 
and the system of equns. 16 and 17 become unstable for 



53 


large values of the weighting matrices. Consequently 
whatever knowledge we already have about the system were 
used^m choosing the initial conditions, it was also 

to give weight age to measurement errors so that 
Ay -no m the behaviour of the parameter estimates in 

approaching the results of the quasilinearization algorithm 
can be used a s a cross check. 


1. The initial condition of the states was 
taken to bo the initial measurement itself, in common- 
with previous methods. 


2. -Starting values for S were chosen about 25 % ■ 
to 30 % higher than the known values, as given by chapter IV. 

3. Taking into account the fact that the matrix 
P is analogous to the error covariance of the states, the 
initial condition P(0) was chosen to be the variance of the 
states and parameters. 


4. The elements of the weighting matrix 0 
(chosen to be diagonal) are determined by systematic trial 
and error, the idea being that for a' chosen initial conditions 
the Parameter estimates should t.end to the known values 
as rapidly as possible and the estimate of the states should 
track the measurements more or less closely. It is found 
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atnx with toe following diagonal elements in 
the given order satisfies the requirements reasonably 

1,611 : 80- °’ 200 -°’ 160 -°- 160 - 0 , 100 . 0 , 100 . 0 , 40 . 0 , 
40 - 0 , 40 . 0 , 35 . 0 , 25 . 0 , 25 . 0 . 

Conclusions 


integration is performed using second 
order predictor-corrector method. The results are shown 
In tables 5.1 and S . 2 . Fig.S.l gives the sequential 
estimation of the parameters 3j_, S 2 , S 3 , 3^ and Spg. It 
Is seen by comparison with table 2.1 that the algorithm 
tracks the measurements closely except near the start. 


To check the reasonableness of the final set of selection 
functions furnished by the method, the program was rerun 
feeding back thefinal set of parameters as initial setting, 
(with corresponding modification of initial variances). 

The weighting matrix was unaltered. This final set of 
parameters is shown in the last column marked RUh II m 
table 5.2. And the estimate for the parameters nearly 
converges to that given by chapter 17. Though the result 
seems to be satisfactory, it has to be conceded that we 
had a good knowledge of the parameters from the very 
beginning which decidedly influenced the initial condition 
settings and weighting factors. 
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To summarise, it is doubtful whether the 
invariant imbedding approach will be successful in a 
similar problem where data is limited, no noise statistics 
are available and knowledge of the parameter values is poor. 
Basically, lack of knowledge about the initial condition 
and weighting factors for the coupled non-linear differential 
equations to be solved is a distinctive disadvantage for 
numerical solution of the same. Also numerical experimen- 
tation with different initial conditions is unattractive 
for a higher dimensional problem. Finally, in the case 
of a system with no driving noise, the basic idea of intro- 
ducing one, and then optimizing an objective function using 
the Minimum Principle does not seem justifiable. 
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Final value of Performance Index = 0.15346 




CHAPTER "fVT 


MAXIMUM LIKELIHOOD IDENTIFICATION 
6.1 Introduction 

One of the most widely accepted methods of 
estimation of non-random parameters in a dynamical system 
under noisy environment is the maximum likelihood esti- 
mation scheme. Indeed, it can be interpreted even in the case 
of random parameters as a form of Baye's estimator where 
no apriori knowledge of the parameter statistics is available. 
The basic idea is that any reasonable estimate of a parameter 
is that value which causes the given observation most likely , 

Ie ? that parameter value which maximises the conditional 
probability density induced on the observations. This 
estimate is called the maximum likelihood estimate. 

Briefly the approach consists of obtaining the 
joint conditional density of the observations given the paramete 
© in accordance with the dynamical model and noise structure 
chosen. 

This gives the likelihood function which is 
subsequently, maximised subject to the system equations. 

When successful, the method yields second order statistics 
of the disturbances also. 

Although in the following derivation it is 
assumed that the noise distribution is Gaussian, it is, 
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rather inappropriate in the sense that in our problem the 
probability of obtaining any negative reading is zero and 
hence the spread of the disturbance has to be restricted. 

But mathematical convenience dictates the choice of a 
Normal probability density function for the noise sequences. 

6.2 Reformulation and Solution of the Problem 

The approach below treats the identification 
problem using the discrete model of the plant and associated 
disturbances. Here we use the matrix A as the discrete 
state transition matrix for a small but finite duration T 
of the process. As in the previous chapters, the process 
is characterised by a set of linear stationary first order 
simultaneous differential equations. It is also assumed that 
the noise sequences are uncorrelated Gaussian random processes. 

Accordingly the. system equations are given by 

x(k) = Ax(k-l) + w(k-l) 

y(k) = x (k) + v(k) k = 2, N (6.2.1) 

where A is the unknown plant matrix to be identified, x(k) 
is the state and y(k) the measurement vectors (both of 
dimension r), v(k) and w(k) are stationary zero mean uncor- 
r elated noise sequences of dimension r. 
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The second order statistics are given by 

E j w(i) w T (j)j = F&J. (6,2.2 ) 

and E I v(i) v T (j)l = US : , (6.2.3) 

l J J 

where 6 -y is the Kronecker Delta, 

=1 when i = j 
tJij 0 when i ji 3 

Also the driving and measurement noise sequences 
are mutually uncorrelated. 

ie, E | v(i) w T (;j)] = 0 for all i, j 

and E [v(i) y T (j)| = 0 for all i, o' 

The Likelihood F u nction " o ■ ' . 

We wish to determine the maximum likelihood 
estimate of the matrix A and the noise variances F and R. 

We denote the elements of A and the elements of B (a r x r 
matrix, which is introduced later) by the vector 9. The 
conditional probability density function 

p fy(k)/y(j), 1 4 3 4 k -1 5 ©] of the measurement at the 

k th stage given all previous measurements and parameters 
is Gaussian by the assumption on the noise statistics. 
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and is given by N fy(k/©), P(k)] where 

y(k/©) = E |y(k)/y(j) 1 5 3 $ k-1 ; ©j 

[ (yOO - y(k/©)}{y(k)-y(k/9)j 


and P(k) = E 


t] 


(6.2.4) 

(6.2.5) 


It may be noted that as more and more obser- 
vations are processed, y(k) tends to the true value and the 

error convariance matrix P(k)' tends to a constant matrix P 

/ 

which is positive definite. 


From the system equations 1 and 2, we can obtain 
a recursive relation for y(k/9) as follows .(Kashyep, 1970) 

3t(k) = A jV(k-l) - v(k-l) j + w(k-l) 
y(k) = A y(k-l) -Av(k-l) + w(k-l)+v(k) 

= Ay (k-1) + w x (k) (6.2.6) 

where wp(k) consists of two noise sequences w and v. 

We define a new noise sequence Wg(k) given by 

W2(k) = e(k) - Be (k-1) (6.2.7) 

where e(k) is the error involved in the estimate of y(k) 
and is a r dimensional vector 

ie, e(k) = y(k) - y(k/©) 
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where E (e(k)j = 0 and E {e(k) e T (j)}= pA (6.2.8) 

Then replacing w 1 (k) by w 2 (k) in eqn. 6, 

y(k) - Ay(k-l) = e(k) - Be(k-l) (6.2.9) 

The noise sequences w-^(k) and Wg (k) have to 
satisfy equality of second moments 

ie E | w^(k) w-j^(k)| = E / wg(k) Wg^(k)^ 

ie ABA T + P + R = P + BPB T 
Therefore, F = P - R + BPB T r- ARA T (6.2.10) 

Also E ^w-j_(k)wj_^(k~l)j- - E / wg (k)wg^(k~l)| 

ie, AS = BP ; Therefore, R = A“^BP (6.2.11) 

We shall subsequently estimate the '‘best” value 
of P and use the same for obtaining R and F. 

The likelihood function for large N can be 
now written as 

% (o, p) = in p [y(D, y(N)/e ] (6.2.12) 

N 

= In p fy (1/9)1 + I In p fy(j)/y(k) Uk ^ h-1 5 ol 

L 3=2 J 

= - rH/2 In - N/2 In | det P | - i f(||y( j)-y(j/©>l 
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rN m 

— In 2TT - _J1 
2 2 


In 



N 

I||e(J/9) 

3-1 


HPp 


-1 


(6.2.13) 


where pfy (1/9)]J is chosen arbitrarily. 

At this stage the problem reduces to one 
of maximising Ljj( 9, P) with respect to 9 where e(k/9) 
obeys eqn.9. 

Estimation of P, A and B 

Choosing nominal starting values for A and 
B matrices we can determine the error sequence e(k), 

k = 2, N. Consistant with the assumption that 

there is negligible error in the measurement of initial 
conditions we set e(l) = 0 to initiate the sequence e(k). 

Maximising the likelihood function Ljj(©,P) 
w.r.t. P we have 


dL N (9,P) 
3 P 


= 0 
P=P(N) 


Simplifying we get 
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P(N) = Best estimate of P using N measurements 
1 — 

= > e(.l/9) e T (j/e) (6,2.14) 

Now substituting back in the likelihood function 

we get 

%(©, P (N) ) = jjL+ln 2T[J In det [p(N)J (6.2.15) 

2 2 

Now maximising R.H.S. of eqn.15 is equivalent to 
minimising det P(N) which is redefined as Jjj(9) 

Hence J N (0) = det [P(N)J 

N 

= det Y e(j/©) e T (j/©) (6.2.16) 

J=1 

Now the maximum likelihood estimates of A and 
B are obtained by minimising Jf[(©) w.r.t. A ana B where 
e(k/9) obeys the eqn» 

e (k/©) - B e(k-l/©) = y(k) - A y(k-l) k=2, ... N 

In addition we have the physical constraint 
on A which must be lower traiangular. 

ie ? a^j — 0 


j > i, i = 1, 2 


r (6.2.17) 
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The procedure now entails evaluation of the 
constrained partial derivatives S Jjj (©)/ d A and £ J N (9)/<£ B 
and is discussed subsequently. 


At this stage a practical difficulty in 
implementing the scheme as above has to be considered.. As 
can be seen from eqn. 16, the matrix P(N) will be singular 
for any N <r. One necessary (but not sufficient) . condi- 
tion for P(N) not to be ill conditioned is that N is very 
large compared to r. This is not satisfied in our problem 
and serious numerical errors crop up in the evaluation of 
det P(N) and ^P(N)’]”^-. Both have to be evaluated repeatedly 
in the iterative procedure. After considerable numerical 
experimentation using double precision, it was concluded 
that for a system of moderate dimension as this, the solution, 
if at all possible will be of little value because of 
computational errors. Accordingly the objective function to 
be optimized was modified to 


N 


J<0) = eT( 3/£) e(j/9) 


3=2 


( 6 . 2 . 18 ) 


which is again minimised subject to eqns. 9 and 17 
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It is not claimed, that minimising the sum 
of squares of error is the same analytically as minimising 
the determinant of error conviance matrix. But under the 
specific conditions imposed by the problem it seems best 
to go in for the alternative, since, intuitively the 
modified criterion seems to be reasonable. 


The augmented function Ja with A(k) and 
as Legrange Multipliers is formed. 

N TvT T 

^ a (©) = ^& T (k)e(k) +52 X ^ [e(k) - B e(k-l) 
k=2 k=2 

- y(k) + A y(k-l)] + a (6.2.19) 

where a is a vector consisting of all super diagonal 
elements of A which are constrained to be zero. 

The A(k) sequence and are evaluated by 

setting -A£a__ = 0 C6.2.20) 

d e(k) 

and — — ■ — — " — — 0 (6.2.21) 

d a 

The first equation gives 
X(k) = -2 e(k) + B^X(k + 1) (6.2.22) 
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Also . ACN+1) = 0 

Hence by integrating backwards A(k), 
k = 2, H are evaluated. 

The gradients 3J a /3A and are 

evaluated by partial diff erantiation taking into account 
eqn 21 also alongwith. 


Therefore, 


N t 

— — - = 27 A(k) y 1 (k - 1) 

$ A k=2 


with 


J a 

■_ d A 


ij 


0 


3 > i 

i = 1, 


(6.2.23) 

r 


and 



j>A(k) e T (k - 1) 
k-2 


(6.2.24) 


This completes the derivation of the essential 
equations for the optimization procedure. 


6.3 Computational Aspects 

Because of the very nature of the formulation 
and solution as discussed above, it is apparent that meaning- 
ful results can be expected only when large data is available. 
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Moreover when measurements are available at non uniform 
intervals the computational algorithm becomes very compli- 
catedy with no closed form expressions available for the 
various gradients. In that case it becomes necessary to 
evaluate them numerically which takes a prohibitively long 
time. So in common with the invariant imbedding approach, it 
was decided to interpolate the data corresponding totfs^. cy 
every 5 revolutions of the mill. Different interpolating 
schemes were considered and were compared with a simulated 
solution using the results of chapter IV. It was observed 
that linear interpolation involves minimum error followed by 
a quadratic polynomial, and hence the former was adopted. 

The minimisation was done iteratively, 
starting with a nominal value for A and B. One iteration 
consists of evaluation of the objective function, the 
Lagrange multipliers A(k), and the gradients w.r.t. A and 
B. An one dimensional minimization was carried out along 
the negative gradient direction by a recursive interval 
halving technique, and the optimum step size determined. 

The starting values are now updated and the procedure 
is repeated. Shopping criterion is less than 2$ improvement 
in the objective function or sum of the absolute values 
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of the elements of the gradient w.r.t. A and B be less 
than 0.30- and 0.50 respectively. 

The program was run with close starting values 
for A. Since B is unknown, it is set equal to A initially. 
The selection functions at the end of 28 iterations are 
shown in table 6.1, It is seen that the convergence as 
far as the selection functions are concerned is quite 
slow. This manifests itself very adversely (in that it 
takes an extravagant amount of computation)when the starting 
guesses are off by about 20% or more. It appears that the 
lack of knowledge of the initial settings of B is mainly 
responsible for this situation. 

One aspect of the method is worth pointing out. 
Since we are identifying the whole A matrix, by determining 
the corresponding continuous system matrix, we can evaluate 
all the breakage parameters without the assumption of 
normalizability . The entire M Z matrix 11 thus evaluated is 
shown in table 6.2 and indicates that the assumption about 
the normalizability of breakage parameter is justifiable. 
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Fig. 6.1 gives the minimisation of perfor- 
mance index J and the monotonic behaviour may be noted. 

The covariance matrices P, R and F are next 
calculated using equations 6.2.10 and 6.2.11 and are shown 
in table 6.3, 6.4, 6.5. 

6.4 Discussion 

From the l, most likely" estimates of the 
covariance matrices P, R and F, it is seen that the driving 
noise is much small compared with the measurement noise* 

Also, the variance of error in the estimated output and the 
measurement error covariance are nearly equal. This indicates 
that the model can be simplified considerably without serious 
error by ignoring the driving noise. It may be noted that 
this is in conformity with the earlier model in chapter III 
and IV which on the whole yield comparable results. 

Due to the slow convergence of the algorithm, it has 
not been found possible to verify whether we will get consis- 
tent results for large errors in the initial guesses. The 
results at the end of 26 iterations for such a case is shown 
in table 6,1. All the same it Is noted that the parameters 
approach the previous values. 
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CHAPTER VII 


STOCK ASTIG APPRQXIMATT f)¥ 

7,1 Introduction 

An interesting approach to the problem of 
parameter identification characterised by a large amount 
of 'dirty 1 data and the absence of any knowledge of 
statistical information is through the methods of 
Stochastic Approximation. The Stochastic Approximation 
algorithms •ften referred to as learning procedures may 
be considered as the stochastic analogs of their determi- 
nistic counterparts. An identification problem in the 

framework of stochastic approximation may be stated as 
follows i 

Let © be the parameter we are trying to 
identify and J(Q) the function to be minimised, which 
however is unknown. But we have experimental results 
where we observe y( 9) and let g(Q) be a known function 
Of y(9) and Q such that E [g (©)/©] = J(9). So 
essentially we look for an algorithm by which we can 
determine © = 9° which minimises J(9) by requiring only 
measurements of the random variable y(©). 

7,2 Ib e Kiefer Wolfowitz Algorithm 

As seen in the last section we are interested 
in finding out the extremum of a function J(©)=E [g (©)/©] 
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where we are able to determine- g(9) given any nominal value 
of the parameter ©. But no knowledge of the statistics 
of the random variable g(9) is available. The Kiefer- 
Wolfowitz scheme uses the estimate of the gradient 
<?©(g (©)/©] computed from neghbouring values of g(©) by 
finite perturbation of © . When © is a vector (say of 
dimension n) one component is perturbed at a time. We 
observe the random samples of g (© ) which are gC^+c^e-) 
j - 1,2 ...... n where c^ is a suitably chosen positive 

sequence, and e ± , e 2 e n are the n dimensional unit 


Vectors : 

e l = 

H 

O 

N* 

• 

• 

• 

* 

... 0 ) 


e 2 = 

( 0 , 1, ... 

... 0 ) 


e n= (°i 0, 1) T 

The subscript ^ indicates a sequence for the 
recursive updating of © as indicated below: 


n 


Thus V© gfetc) 


SI e - ~_ 

j=l ° 2 ^k 


g (©k+cjj-e-j ) -g (©j^-c^ej ) | 


(7.2.1) 
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Using the estimate of the gradient, the 
computing scheme can be written as 

G k+l e k - V 9 £g(©k)] ( 7 . 2 . 2 ) 

where d k is a suitably chosen sequence. 


For convergence, the sequences c k and d k 
have to satisfy the following conditions: 


^ d k ~ 


OO 


£ a k ° k ' °° 


Lt Cjj = o 

J£ “7* 00 


d k c di c k } 0 < d-L < 


OO 


and (d k / c k ) 1 2 3 < c*o 


(7.2.3) 


The conditions that should be satisfied for 
convergence of KW scheme as given by Venter (1967) are: 

1. J(9) and its second order derivatives are 
bounded in R n , 

2. ©° is a local minimum of J(9). i e , for 
some £ > 0, we have 

J(9°)^J(9) for all © = 9° in the set |© : j{ © . ©o j| c 

3. For every © € R n , V Q J( 9 ) j£ 0 if © ? ©? 
That is ,© 0 is the only stationary point of 1 ( 9 ). 
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4. When J (9) is redifined a s J (©) = g (Q) + v (9) 
with E{v(9)l = 0, v(9) should satisfy the condition 

B{v 2 (9)}<-o ¥ © £ R n . (7.2.4) 

Venter has shown that if these conditions 
are met, the algorithm given by (7.2.2) with gains satis- 
fying (7.2.3) converges with probability one to either 
or infinity. 

7.3 Computational Considerations 

In applying the algorithm to our problem, 
g(©) is chosen as the sum of squares of error between 
the nominal trajectory given by the parameter value 9 
and the experimental data. This enables comparison with 
the results of the previous chapters. 

The sequences cq- and dq- are chosen so as to 
satisfy (7.2.3). The usual choice is dk = dq / k and 
°k = Ci/k^where dq and Cq are positive and 0 < ^ < 0,5_. 

We choose < = 0.4, and cq adjusted so that the perturbance 
on the performance index is neither too largenor too small, 
dq was chosen by trial and error to get rapid convergence, 
and range from 0.5 to 40.0 for the 12 parameters. It is 
seen that apart from the rate of convergence, the choice of 
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d x within reasonable limits does not affect the eventual 

result seriously, provided there is no violent oscillation 
for any parameter values. 

It is noted that although the algorithm is 
quite simple, it is computationally inefficient. For each 
perturbation of every parameter, the nominal trajectory 
has to be evaluated and change in performance index 
computed. Hence it was decided to identify only the 
selection functions by this method, keeping the breakage 

parameters at the value given by the Quasilinearization 
method. 

7 *4 Results 

Table 7.1 summarizes the numerical results 
of the estimation of the selection functions for two 
different starting values. It is seen that by adjusting 
the coefficients d^, It is possible to obtain faster 
convergence even if the initial settings are far off 
from the optimum values. The minimization of the performance 
index is also given. It can be observed that the results 
for thu, two different starting guesses are more or less 
consistent. The convergence of and $2 are shown in 
figure 7.1. 
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CHAPTER VIII 


CONCLUSION 

The problem of identification of a dynamical 
system in noisy environment is essentially that of stoch- 
astic observability of a set of unknown parameters, given 
a set of related measurements. Although several procedures 
have been suggested in the literature for the identification 
of systems, in the absence of any explicit creiteria for 
the observability of the parameter, the feasibility of 
these approaches must be judged in the light of their 
computational complexity, the amount of numerical experi- 
mentation required and the trend towards convergence of 
the estimates for a given data. Clearly any conclusion 
regarding these factors has to be restricted to a specific 
problem. It is in this premise that this thesis is 
investigated where an attempt has been made to study the question 
of identification of a batch grinding mill - a problem of 
considerable engineering importance. 

Different approaches for the identification 
of the given system were discussed in the preceeding 
chapters. A comparitive evaluation of the procedures is 
considered appropriate. At the same time the question as 
to whether comparisons will be meaningful has also to be 
examined, since every method has one or other distinctive 
aspect associated with it. 



83 


As far as the performance index is concerned, 
there is close similarity between the methods of Recursive 
Lease Squares, Quasilinearization and Stochastic Approxi- 
mation. So it is natural to expect comparable results 
from the three schemes. As far as the selection function 
goes, it is seen that this expectation is satisfied to a 
large extent. It is also worth pointing out that in all 
the three methods, the initial guesses for the parameters 
is not a serious issue. Except that more computer time 
will be involved, a large error in the initial setting 
does not materially alter the final results. A still more 
important observation relevant here is that the available 
data which is limited and of non-uniform distribution 
does not pose any major difficulty either in the implemen- 
tation or convergence of the program. It was not considered 
necessary to resort to any formof interpolation. As far 
as the specific problem is concerned, it is seen that 
Quasilinearization promises to be most suitable in many 
respects. Identification of both selection and breakage 
functions is achieved without lot of involved computation. 
Convergence is mono tonic and rapid. 

The method of Invariant Imbedding stands apart 
in many respects. Briefly the identification scheme calls 
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for more information about the process than is available. 

It is possible that numerical experimentation can answer 
some of the specific issues like choosing weighting matrices 
etc satisfactorily in a low dimensional problem, but not 
necessarily so otherwise. It is in fact possible to pose 
the question of best choice of weighting matrices in a 
situation where noise statistics are nor available as a 
problem in itself, and has not been answered here. On 
the whole it is felt that the results of this chapter only 
marginally complement rather than supplement what is 
already known. 

The maximum likelihood approach differs from 
the rest in that, we seek to identify the discrete plant 
matrix in full, and estimate the associated noise convariances . 
A more comlete description is thus obtained. The performance 
index is apparently the same as before, but is evaluated in a 
different way. The estimate of the covariance cf the 
driving and measurement noise is considered significant, 
since it enables comparison of the relative strengths. As 
has been mentioned earlier, the assumption of Gaussian 
distribution is ,ina-PP r °priate in our problem. It will be 
interesting to see whether some non Gaussian density can be 
employed, but it is likely to introduce mathematical compli- 
cations . 

li 
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Table 8.1 summarises the important results 
of the different schemes for easy comparison. 

Briefly the main results of this thesis are 
the following: 

1. The linear time invariant model that was 
chosen appears to be quite adequate to describe the system 
behaviour in the time interval of interest f 

I 

2. Since no general results concerning 
convergence or uniqueness of parameter values are obtained, 
the fact that we are able to obtain reasonably close 
estimates using different procedures is significant. 
Heuristically speaking this solves the problem more or less 
satisfactorily. 

It might be possible to increase the accuracy 
of the computations by using more sophisticated techniques 
like higher order predictor-corrector schemes for solving 
the differential equations, improved methods for matrix 
inversion etc. Also approaches like instrumental variables 
or time series analysis can considered. But basically it 
seems there can be no escape from the same sort of difficulties 
as above, posed by moderately large number of parameters to 
be identified with limited data coupled with lack of perti- 
nent statistical information. 
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